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Abstract 

We derive a new formulation to calculate the excess chemical potential of a fraction of A'^i 



> 

(N 

00 

f— N I particles interacting with A'^2 particles of a different species. The excess chemical potential is 

\0 ■ 

C^ ' calculated numerically from first principles by coupling molecular dynamics and Thomas-Fermi 

^' 

^D ■ density functional theory to take into account the contribution arising from the quantum electrons 



on the forces acting on the ions. The choice of this simple functional is motivated by the fact 



O I that the present paper is devoted to the derivation and the validation of the method but more 

^ ■ complicated functionals can and will be implemented in the future. This new method is applied 

in the microcanonical ensemble, the most natural ensemble for molecular dynamics simulations. 
This avoids the introduction of a thermostat in the simulation, and thus uncontroled modifications 
of the trajectories calculated from the forces between particles. The calculations are conducted 
for three values of the input thermodynamic quantities, energy and density, and for different total 
numbers of particles in order to examine the uncertainties due to finite size effects. This method 
and these calculations lie the basic foundation to study the thermodynamic stability of dense 
mixtures, without any a priori assumption on the degree of ionization of the different species. 
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I. INTRODUCTION 

The thermodynamic stabihty of dense ionic mixtures bears important consequences not 
only on our understanding of the thermodynamic properties of dense binary systems but also 
on the structure and the evolution of gaseous giant planets. Indeed, the interior of jovian 
planets is composed essentially of hydrogen and helium under either atomic or molecular 
form in the outermost envelope and under the form of a partially or fully ionized plasma 
in the inner regions [l| j^. Temperatures and pressures along the Jupiter or Saturn in- 
ternal isentrope conditions range from about 100 K to ~ 20000 K and from about 1 bar 
to ~ 60 Mbar. Under these conditions, not only the hydrogen/helium mixture experiences 
pressure ionization but the homogeneous mixture may become thermodynamically unstable. 
Such an immiscibility between helium-rich droplets and a hydrogen-rich fluid will liberate 



extra gravitational 



energy, modifying significantly the energy balance and thus the cooling 



of the planet J3| J4| [5|. For terrestrial applications, inertial confinement fusion experiments 
or laser-driven shock-wave experiments on hydrogen isotopes reach densities and tempera- 
tures characteristic of the aforementioned planetary conditions, probing the thermodynamic 
properties of dense plasmas and requiring a correct theoretical foundation to describe their 
equilibrium properties. 

The theoretical description of the thermodynamic phase diagram of a dense two- 
component system is a particularly complicated task, for it requires a correct description of 
the excess free enthalpy (in a pressure-temperature diagram) of the mixture with respect to 
the pure phases. This excess quantity is very small compared with the contributions of both 
the mixture and the pure phases (it is by definition close to zero near the critical point) and 
therefore must be calculated with very high accuracy. Early calculations, based on simplified 



analytic or semi-analytic calculations of the 



ree energy of the plasma, assumed hydrogen 



and helium atoms to be fully ionized p 7| 8[ J9|. Moreover, these calculations assumed 
either a rigid electron background, the so-called binary ionic mixture (BIM) model, or a po- 
larizable electron background within the linear response approximation. Although correct at 
very high density or temperature, these assumptions fail when electrons and protons start to 
recombine. The phase diagrams calculated under these conditions are thus restricted to a re- 
duced (high) density-temperature range. Further attempts to do a consistent, first-principle 
determination of the H/He phase diagram, with no assumption on the electron distribution 



around the ionic centers and a correct treatment of the various A^-body ion and electron 
interactions, were based either on extrapolation at finite-temperature of zero-temperature 



calculations [10'] or on incorrect thermodynamics integration jll| and thus remain also of 
doubtful validity. Under such circumstances, it is clear that not only the thermodynamic 
phase diagram of a hydrogen-helium system at high density has not been established accu- 
rately yet, but the correct calculation of the excess free enthalpy of a concentration of atoms 
immersed in an interacting system of different species remains to be done. 

In this paper, we derive a new method to address this very point, which is crucial for a 
reliable determination of the thermodynamic phase diagram of dense binary mixtures. This 
method is applied in the microcanonical ensemble and allows the direct calculation, from 
first principles, of the excess chemical potential of a binary mixture of nuclei and electrons 
interacting through the Coulomb potential. Calculations in the microcanonical ensemble 
allow a fully consistent calculation between the forces acting on the particles and the induced 
trajectories, without the introduction of thermostats. We first derive the thermodynamic 
equations which allow the exact determination of the excess chemical potential. We then 
combine the density functional theory (DFT) to describe the quantum mechanical properties 
of the electrons and molecular dynamics (MD) to integrate the ion classical equations of 
motion to calculate this chemical potential. Since the present paper is devoted to the 
derivation of the method, we use a simplified functional form for the electrons, namely 
the Thomas-Fermi approximation, in order to speed up the minimization of the energy. 
We limit our calculations to three different values of the appropriate input thermodynamic 
quantities, namely energy and density in the microcanonical formulation used in the present 
paper. In a future work, devoted to the global analysis of the H/He mixture under various 
thermodynamic conditions, a more general functional form will be implemented. Section 2 
presents the derivation of the chemical potential of A'^i atoms of a given species interacting 
with N2 nuclei of a different species. Section 3 describes our general energy functional 
to take into account the quantum behaviour of the electrons when computing the ionic 
configurations, a necessary condition for an accurate treatment of the problem. Section 4 
is devoted to the description of the MD numerical computations, to the discussion of the 
finite size effects and to the presentation of the results obtained for different thermodynamic 
conditions. The last section is devoted to the conclusion. 



II. DERIVATION OF THE EXCESS CHEMICAL POTENTIAL 

As mentioned in the introduction, the ultimate goal of our calculations is to determine the 
thermodynamic stability of a given number A^^i of atoms immersed in a system of N2 particles 
of a different species under given thermodynamic conditions, without any assumption on the 
electron distribution around the nuclei, i. e. on the degree of ionization of the atoms. The 
stability of such a mixture involves the calculation of the mixing enthalpy of the system, 
i.e. of the excess chemical potential of each immersed atom. The chemical potential /ij of a 
particle i immersed in a plasma corresponds by definition to the change of the state function 
of the appropriate thermodynamic ensemble when one adds or removes this particle to/from 
the plasma. When the thermodynamic limit is achieved (A^ -^ 00, V -^ 00, N/V =constant), 
the result does not depend either on the ensemble or on the fact that the particle has been 
added to or removed from the surrounding plasma. Because of the large fluctuations of the 
system away from its equilibrium configuration, one can not add or remove directly a particle, 
in particular at high density or if the interaction potential is too stiff. The correct approach 
consists in modifying progressively the interaction potential A V{r) between the particle i and 
the surrounding particles j 7^ i. The case A = corresponds to the case where the particle i 
does not interact with the other particles, but retains its discernability character (case of an 
ideal mixture), whereas the case A = 1 corresponds to the sought two-component system with 
full interactions. This method illustrates the so-called thermodynamic integration approach. 
In the microcanonical ensemble, with fixed energy, volume and number of particles {E, V, iV), 
the chemical potential of a particle "1" of mass mi corresponds to the calculation of the 
following expression derived in appendix A-C: 
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E{X)-V d\ /' 

where F is the Gamma function, V is the cell volume, Kj<q is the kinetic energy and E{\) is 
the energy corresponding to the system with the interaction potential A V{r). The brackets 
(...) denote a microcanonical average. The first term on the right hand side is the ideal 
part of the chemical potential, arising from the entropy cost due to the particle insertion 



or removal, while the second term represents the non-ideal contribution of the chemical 
potential, which depends on the interaction between the particle under consideration and 



the rest of the system. The integral can be estimated by a Gauss-Legendre quadrature 
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where Xi are the zeros of Legendre polynomials and tOi are the associated weights |l^ . 

These calculations, in practice, require great caution. By definition of the microcanonical 
ensemble, the total (potential+kinetic) energy of the system must be conserved along the 
simulation. As a consequence, at each step where the potential goes from \iV to Aj+iV^, 
the kinetic part of the energy must be renormalized in order to maintain the total energy 
constant. Therefore, the correct calculation of the chemical potential consists in generat- 
ing several particle configurations (to be described in §4) corresponding to the Hamiltonian 
H + XV, and in computing the averages which appear in Eq.(0). This can be done by a 
fully classical simulation if the interaction between each kind of particles is described by 
a classical 2-body potential between particles. Such an approach, however, can not take 
into account the fact that the potentials strongly depend on the density and the tempera- 
ture, evolving from a potential characteristic of atoms at low density and temperature to a 
long-range Coulomb potential characteristic of a fully ionized plasma at high density and/or 
temperature. Therefore, the correct phase diagram, without any assumption on the inter- 
action potentials, requires ab initio generations of representative ionic configurations. This 
approach is presented in the next sections. 

III. THE FUNCTIONAL OF THE ELECTRONS 

As mentioned in the introduction, a correct study of the problem under consideration 
requires a correct treatment of the ion and electron interactions. This implies to take 
into account the effects of the quantum nature of the electrons on the forces acting on 
the ions. Since the pioneering work of Hohenberg and Kohn [14], many problems involving 
interacting electrons have been tackled within the framework of the density functional theory 
(DFT). This theory turns the diagonalization problem of a many-body Hamiltonian into 
the minimization of a functional fi[n(r)] of the electron density, a much easier approach 



when dealing with many electrons. For this reason, the DFT has been extensively used 
in condensed matter and is described in detail in many textbooks (see e.g. [l5|). In the 
framework of the DFT, the grand potential of the electrons can be written in the form: 

n[n{r)]= f dr (VUr) - fi)n{r) + F[n{r)l (3) 

where F[n{r)] is a universal functional of the ground state density n{r) of the interacting 
electrons and Vie{r) denotes the external ion-electron potential. fi[r7,(r)] is minimum when 
n{r) corresponds to the correct density. In our calculations, we have chosen to write i7[n(r)] 
under the following simplified form (in order to speed up the minimization): 

Q[nir)] = ksT [ drn{r)^^ + c,^ [ dr[n{r)f' + E,,,,[n{r)] 

f T Nn(r)f e^ f f ^ , n(ri)n(r2) 
+ cw / dr^ — ^-^ + ^ / / dridrs- ^ ^^ ^ ^^ 



n{r) '2 J J 1^1 ~ ''"2I 

drn{r) {VUr, Ri^n) - yu) , (4) 

where Cgx = — |(3/7r)-'^/^e^ is the exchange Dirac coefficient, cw = icr/8)h^/me is the von 
Weizsacker gradient correction coefficient [l^ [la] (with a = 1 in our case), and -E'corr['^(^)] 
is given by a parametrization [l2| of Monte Carlo simulations jlS|. -F3/2(^) and Fi/2{r]) are 
the Fermi integrals, where 77 = (/x — V{r))/kBT is obtained by the inversion of the relation: 
n{r) = (2n'^)^^(2me/h^Y^'^{kBT)^^'^Fi/2{r]). Accurate fitting formulae of the Fermi integrals 
and inverse integrals have been published in the literature |l9| J20|. For a given configuration 
of the nuclei, we are able to find the electronic density n(V) which corresponds to the ground 
state of the system. The Hellmann-Feynman theorem |2l| J2^ enables us to calculate the 
forces arising from this electron density acting on the nuclei and the stress tensor on a cell 
(of which diagonal terms correspond to the pressure components). 

Within the Born-Oppenheimer approximation, we can thus make a classical molecular dy- 
namics (MD) simulation of the nuclei sub-system, while taking into account in the calculation 
of the forces the quantum behaviour of the electrons. The Born-Oppenheimer approximation 
expresses the fact that the electrons respond instantaneously to a change of configuration 
of the ions, a fairly good assumption for dense ionic systems. The calculation of the last 
term on the right hand side of Eq.(0]), i.e. the interaction with the external ionic potential, 
involving effective pseudopotentials, is described below. 



IV. MOLECULAR DYNAMICS 



Method 



We have computed Eq.(0) for a number of A^i helium atoms of nuclear charge Zi = 2 
and mass M\ immersed in a system of A^2 hydrogen particles of charge Z2 = 1 and mass M2. 
The thermodynamic averages in Eq.(^ are estimated by generating a set of representative 
configurations of the system in a cubic reference cell of size L with periodic boundary 
conditions. This is done by a dynamical simulation of the equations of motion for the ions: 

^^^ = ^. (5) 

where Mj is the mass of the ith nuclei. The forces Fi between particles (or equivalently the 
total potential Viv)) arising from electron and ion A^-body interactions, beyond any linear 
approximation for the electron-induced screening effects of the core potential, are calculated 
from a density functional approach. These forces involve the ones arising from the quantum 
electron distribution obtained from Eq.(j3]) and the ones derived from the interionic potential 
ZjZnC^ l\Ri — Rj\. The equations of motion are solved with a standard Verlet velocity algo- 
rithm [2^. The crucial point of the present paper is that these calculations are completed 



in the microcanonical ensemble, i.e. at constant energy, volume and total momentum J24|. 
Standard simulations, in other thermodynamic ensembles, imply the introduction of a ther- 
mostat, either by reinitializing the velocities "periodically" or by introducing new degrees 
of freedom. These thermostats, however, yield a perturbation of the trajectories, which no 
longer represent the ones determined by the forces. Such unphysical effects are avoided in 
the present microcanonical calculations, which insure full consistency between the forces and 
the trajectories. 

The ah initio calculations, with the aforedescribed functional, have been performed with 
the ABINIT code J25[ . We replace the bare Coulomb potential of the nucleus by a pseudopo- 
tential, which differs from the true Coulomb potential below a cutoff radius rioc, removing 
the cusp constraint at r — i> and avoiding the 1/r singularitv. The pseudopotentials used in 
our simulations are those of Hartwigsen et al. for helium ^^ , and Goedecker et al. for hy- 



drogen 



27[. These pseudopotentials are constructed so as to reproduce with high accuracy 



the Kohn-Sham free energy. 



The aforementioned cutoff radius determines an upper bound in density for tlie domain 



of validity of the pseudopotentials. The ones used in the present calculations 



have 



Hoc = 0.2 bohr, which implies a density limit: a>Tioc5 ^-c rs>0.2, where Vg = a/a^ is the 
density parameter, a = { ^J^/^ Y^^ is the mean distance between nuclei and Oq is the Bohr 
radius. This condition corresponds to p<335 gcm~^. 

As mentioned earlier, in order to maintain the total energy constant during the process of 
switching on or off the interaction, the kinetic contribution E]^^^ must be renormalized. For 
the thermodynamic conditions of our runs, this corresponds to a decrease of i^km because 
the potential energy increases as the interaction is switched on (A = ^ 1). This implies 
a large initial kinetic energy. The condition is more easily fulfilled if one chooses to switch 
off the interaction instead of switching it on (A = 1 — ;> 0). Indeed, during such a process, 
the kinetic part must be increased (instead of decreased), which is always possible. As a 
consequence of this renormalization, we can not associate an accurate temperature until the 
thermodynamic limit is reached. This process is represented on Figure which displays 
the kinetic energy during the whole process, and shows the discontinuities appearing in the 
mean value of i^kin when A changes from Aj to Aj+i. 

B. Results 

We have tested our procedure on a system consisting of 63 hydrogen nuclei and 1 helium 
nucleus. The thermodynamic conditions of our microcanonical simulation are: i^tot = 132.21 
hartree and V = L^ = 57.906 bohr^, which correspond to r^ ~ 0.6, T ~ 2 10^ K and 
P ~ 7 10^ GPa. The reference cell of the simulation assumes periodic boundary conditions, 
with one particle exiting the cell on one side replaced by one entering the opposite side. 
In order for the final results not to depend on the initial distribution, which corresponds 
to a random distribution of the positions and velocities (obtained by a classical molecular 
dynamics simulation to prevent atoms to overlap), we let the system relax during 4000 time 
steps, a very conservative limit for the considered densities and temperatures. Even though 
the main contribution to the total energy at r^ = 0.6 comes from the nearly uniform electron 
background, the forces depend partly on the non-uniformity of this electron density distri- 
bution. Therefore, in order to calculate the forces correctly, the electronic density, more 
precisely the departure of the density from a homogeneous distribution, must be calculated 
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with very high accuracy. In order to fulfil this condition, we require the energy to converge 
within \AE /E\ < 10~^. Unfortunately, high accuracy in the functional minimization does 
not preclude energy fluctuations during the simulation due to the discretization of the New- 
ton equations and, most importantly, to finite-size effects. These points are examined below. 
The equations of motion are solved using the Verlet algorithm with a time step equal to 
dtion = 0.25 a.u. This time step enables us to resolve the dynamics of our system accurately 
for any value XV oi the He-H interaction. Figure El displays the conservation of the total 
energy obtained in our simulation with this time step. The integral Eq.Q is first calculated 
with a M-point Gauss-Legendre quadrature. After the first 4000 time steps to let the system 
relax, an other 4000 time steps simulation is ran to generate several configurations. The 
argument dxE{X)/{E{X) — V) is calculated numerically, by calculating E{X ±0.01) every 10 
time steps for a fixed configuration and E{X) —V. A last run is devoted to the calculations 
of the ideal part of the chemical potential by generating 10000 different configurations of 63 
H and the evaluation of {KJ ). 

We have tested the validity of the M-point quadrature to estimate the integral Eq. (jSJ by 
doing similar calculations, for the same thermodynamic conditions, with a 3-point, 6-point 
and a 9-point quadrature. The results are shown on Figure El and the resulting evaluations 
of the integral are given in Table lU As seen in this table, a 6-point quadrature is enough 
to calculate accurately the integral (j21). Our calculations of the chemical potential — -^ 
of a helium atom embedded in a 63-H plasma for our thermodynamic conditions (— -^ 
corresponds to the ideal part of the chemical potential and — -^ to the excess contribution) 
yield: 



Pi 
kT 


= 14.02, 


kT 


= 5.14, 


kT 


= 19.16. 



(6) 



In order to estimate the A^-dependence of our result, we have also calculated the entropy 
cost which corresponds to the removal of 2 He-particles surrounded by 126 H, and 4 He- 
particles surrounded by 252 H, for the same thermodynamic conditions (density and energy) 
as for the {1 He, 63 H} system. These computations are much more time consuming (the 



computation time scales roughly as t oc N^), and the removal of the helium atom must be 
done 2 and 4 times, respectively. We expect a very small dependence of the results on the 
size (A^) of the simulated system. Indeed, we do not calculate the chemical potential of one 
He atom but the chemical potential of a constant fraction {x^e =1/64) of He in a H-He mix- 
ture. As a consequence, the A^-dependence of the results stems from the interaction between 
a He atom with its replicas (due to the periodic condition boundaries) and from the fact 
that the accessible phase space increases with the total number of atoms. The first effect 
becomes important when the characteristic length of interaction is of the same order as the 
simulated box, i.e. at much higher density than the present simulations. Quantification of 
the second effect requires simulations with different values of A^. The results are presented 
in Table IHI The chemical potential is estimated with a centered scheme, i. e. the chemical 
potential corresponding to a Helium fraction xhc =1/128=0.008 is evaluated from the en- 
tropy difference between the {0 He, 63 H} system and the {1 He, 63 H} system (with full 
interaction). The statistical uncertainties on —^i/kT are ±0.02 for the {1 He, 63 H} and 
{2 He, 126 H} systems, and ±0.04 for the {4 He, 252 H} one (achieving the same statistical 
uncertainties scales as A^^). The results between the three systems for a He fraction equal to 
0.008, as shown in the Table, are thus fully compatible, and no statistically-significant trend 
appears. The same simulations yield also the estimation of the Helium chemical potential 
for different number fractions. All the results are given in Table |Hj and are compatible 
within the statistical uncertainties. For the {1 He, 63 H} mixture, we have also conducted 
calculations for two other thermodynamic conditions, displayed in Table IIIII As expected 
intuitively, it is easier to add an atom in a low density plasma than in a high density one 
(at constant total energy), or in a cold plasma than in a hot one (at constant density). 

It is interesting to compare our results with the limit of rigid electronic background at 
high density for the binary ionic mixture, the so-called BIM limit [29| J30|. Our reference 
conditions, i^tot = 132.21 hartree and V = L^ = 57.906 bohr^, i.e. r^ ~ 0.6, correspond to 
T = 2.2 10^ K, P = 6.7 10^ GPa and AS = 19.15 ^b (Table HH) in the simulation. For this 
density and temperature, the BIM corresponds to P = 7.2 10'^ GPa and AS = 21 fc^. We 
have also ran a simulation at higher density, namely r^ = 0.3, which is close to the density 
limit of our pseudopotentials. The total energy is equal to -Etot = 672.76 hartree. In that 
case, the present calculations yield T = 4.5 10^ K, P = 2.15 10^ GPa, whereas the BIM 
results are T = 4.5 10^ K, P = 2.2 10^ GPa. The small differences between the simulations 
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TABLE I: Integration of Equation |21 



3 points 


6 points 9 


points 


trapeze 


J^^o \e{X)-V dX I 


0.05578 0.05548 


0.05565 


TABLE IL Finite size effects 


on the chemical potential. 




System {1 He, 63 H} 


{2 He, 126 H} 


{4 


He, 252 H} 


XHe = A^Hc/(A^H + A^Hc) 0.004 0.008 0.012 


0.004 0.008 0.012 


0.004 


0.008 0.012 


f^ X 19.16 X 

kT 


19.15 19.15 19.15 


19.17 


19.14 19.21 



and the BIM reflect the contribution due to the electron gas polarization (inhomogeneous 
distribution), which starts playing a role around these densities, and the contribution due 
to the interactions between particles of different species (namely H and He). These latter 
are not taken into account in the BIM, which is based on the so-called linear volume law, 
where only the ideal entropy of mixture is included. 

V. CONCLUSION 

In this paper, we have derived a new method, based on physics flrst principles, to calculate 
the excess potential of a number fraction of particles immersed in a mixture of particles of 
different species, as given by Eq.iH)). The calculations are based on a consistent treatment of 
the forces acting on the nuclei, taking into account the contribution arising from the quantum 
electrons, by calculating self-consistently the equations of motion of the classical nuclei and 
the functional density of the electronic distribution. The method is applied directly in the 
microcanonical ensemble, avoiding the use of a thermostat, and thus insures consistency 

TABLE HI: Chemical potentials for different thermodynamic conditions. 

{£^tot(hartree), ^^(bohr^)} {132.21,57.91} {132.21,139.40} {20.06,139.40} 

r, 0.6 0.8 0.8 



kT 



19.16 16.59 13.94 
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between the forces and the trajectories of the particles. The bare Coulomb potential is 
approximated at short distances by pseudopotentials which remain valid up to large densities 
(rs>0.2), where the linear response theory becomes valid. The thermodynamic quantities 
are calculated for different configurations, representing the evolution of the interaction, 
and thus of the system, from the initial case of an ideal atom "1" inserted in a system 
of particles "2" to the final case where all interactions between the immersed particle and 
the surrounding nuclei are taken into account. Only properly following such a series of 
changes of equilibrium states insures thermodynamic consistency and thus allows a correct 
evaluation of the energy and pressure contributions to the excess chemical potential of the 
immersed particle. Previous simulations [llj calculated the excess enthalpy directly from 
the difference beween the final and the initial states, yielding an incorrect evaluation of the 
contraction work, and thus of the pressure contribution. 

The validity of the method has been tested with the case of a dense hydrogen/helium 
mixture for three different helium fractions and three different thermodynamic states. The 
forces are calculated with very high accuracy, with a convergence criterium \AE/E\ < 10^*. 
Finite size effects on the final results have been quantified and found to be small (~ 10"'^) 
leading to fluctuations of the same order on the total energy (see Fig 2). The method provides 
robust foundations for accurate evaluations of the excess thermodynamic quantities of dense 
binary mixtures, without any assumption on the electron density distribution and thus on 
the degree of ionization of the atoms. This opens the door to accurate calculations of phase 
diagrams of dense mixtures of atoms and partially or fully ionized plasmas, a subject of 
prime interest for the structure and the evolution of giant gaseous planets. Work in this 
direction is in progress. 

We are very grateful to Gerard Massacrier and Alexander Potekhin for very useful dis- 
cussions and insightful remarks. We are also indebted to the two anonymous referees for 
their valuable comments which helped improving the initial manuscript. 

APPENDIX A: MICROCANONICAL AVERAGE 

We consider a (classical) system with fixed total energy E, volume V and total momentum 
Ptot, which contains two different kinds of particles, Ni, N2, with N = Ni + N2. In this 
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microcanonical ensemble, the number of accessible states for this system is: 



6Q 



6E 



h^^N- 



^ /"dp^dg^ S{E - H) 5(p,ot - J2p^), (A1) 

.7=1 



where H = Ylj=i ^m' ~^ ^(.Q^) is the Hamiltonian of the ionic centers, and includes 
the modification of the Coulomb potentials due to the electron gas polarization, dp^ = 
Y[i=i Y[j=idPijy ^iid dpij is the j-component of the momentum of the particle i. Idem for 
dq'^ . Then: 

N 



" = Ssra / ''P"''''" "(^ - "1 *(P« - Eft) - j;w^, m P.O.), (A2) 

i=i 
where d[x) = 1 if a; > 0, otherwise, and J(i?, ptot) denotes the integral. 

The Laplace transform (toward E) of I is J31 1 : 



/p+oa 
dp^dq^ 5(ptot - Yl P^) / d^ ^(^ - H) 

N 

= /dp^dq^<5(p,ot-5^P,)-e-^^ 
-^ i=i ^ 

dp^ 5(ptot - J2pj)^ ""^ J dg"^ -/-' 



j=i 

The Hamiltonian if is general and does not have, in particular, to be positive, although 
it needs to have a lower limit. Note that with a change of variable ( = E — min(if), the 



integral f-^^^/j^-. becomes J^ °°, and the results remain unchanged. 



Then: 



Jo = / dp^ (5(ptot - X^Pj)e 

-^ .7 = 1 



^ N P. 






13 



r? 

N N -sE,=i -^^ 



where J = J dp{^6{ptot^ - Ei=iPiJe '^' ' ""^ • 
J can be calculated by Fourier transformation: 



J^[J] = / dpf — = / dptoi^5(ptoi^ - y^PjJ exp -sVtt^ exp«p, 
I ^ [ ( V^- \ 



m^vA-'^^' 



Then: 



nJLi V27rmj 1 

JV 






"J 



With ptot = 0, we get: 
and: 



3 _ I n7=i V27rmj 



27rEL"^i/ s 



3(iV-l) ' 






Equation (A2) thus reads: 



^^"^i'^2n./2.E,li-.7 -^ r(5(^ + i 



(A3) 



In a similar way, we have for ci;= (|^) 



N2 



For a quantity depending only on q^, z.e. A(g^), we have: 

N 



<-^> " ^Pv7^/dp"d," ^(£ - H) S(p„ -Y,P,} A(q- 



Pj) ^IC^^ 

/ dq^ eiE - V{q'')){E - V{q^))'-^^-^A{q'') 



J dq^ 9{E - V{q^)){E - V{q^)) 
14 



3(iV-l) -I 
2 -^ 



(A5) 



where (...) denotes a microcanonical average. 



APPENDIX B: CHEMICAL POTENTIAL BY THE PARTICLE INSERTION 
METHOD 



The definition of the chemical potential of the particle 1 is: 

/ii / dS 



kT 



dNi 



(Bl) 



E,V,N2 



With S = klnu and the equations derived in Appendix A, this yields: 

_ }M_ _ lnu;Ar^+i -Inc^Af^ _ Un^+i 
kT ~ Ni + l-Ni ~ ^ 



In 



2711711 



^Afi 



'^'^EUrn, 



3/2 



r ( ^^j J^qN+l 0^E _ ^)(^ _ y^^N+i^y-f-i 



ElX'rnJ N, + l r(f) fdqNo^E-V)iE-Viqnf-^-^ 



:b21 



Let define I^ as: 



'Af 



3(iV-l) 



7VxxH^^_l 



dq^e{E-V{q^))K-^ 

where K^ = E — V(q^) is a function of q^ and should not be formally confused with the 
kinetic part of the Hamiltonian (even if Kn is equal to the kinetic energy in a molecular 
dynamics simulation). 
We get: 

dg^+i y"dg^ e{E - V){E - \/(0 - ViqN+i))'^-' 



IN+I 



dq 



N+l 



3 3(iV-l) -, 

dg^ e{E - V)KlK^ ^ II- 



V{q 



N+lf 



K 



N 






and: 



'Af+l 



dg 



Af+l \ ^^N 



Kld-''^^--^ 



K 



N 



which yields for the chemical potential: 



kT 



In 
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3/2 -p ( 3{N-1) 

1 ^ I 2 



^J iVi + 1 r(f) 



dqr^+i ( i^^ I 1 - 



^(g_ 



Af+ly 



is: 



TV 



(B3) 
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APPENDIX C: CHEMICAL POTENTIAL BY THE THERMODYNAMIC INTE- 
GRATION METHOD 



The variation of entropy when going from a state with interaction A = to A = 1 reads: 

Jx=o d\' 
We can thus derive: 



AS= dX — . (CI) 



]_dS_ _ l_du 
k d\ 00 d\ 

_ dxj dq^+i ^(^(A) - \/(q^+i))(E(A) - V{q''+^)y 

~ /dg^+i e{E{X) - V{q^+^)){E{X) - V{q^+^))^~^ 

3N \ /dq^+^ ejEjX) - V-(q^+^))^f (i^(A) - V-(g^+^))- 

2 ; / dg^+i e{E{x) - v{q^+^)){E{x) - viq^+^yi'f- 

3N \ / 1 dE 



2 J \E-V dX 

The thermodynamic integration proceeds in two steps. The first one deals with the insertion 
of a free particle into the system. The entropy cost of this insertion is given by the formula 
established for the insertion method. 



^ = In 

kT 






V h^ J \El\'rnJ N, + l r(f) 



(C2) 



where V is the cell volume. 

The interaction is then progressively switched on, and the non-ideal part of chemical poten- 
tial is then given by: 

kT 



^^c^-^c-im^'m- 



The total chemical potential is the sum of the two contributions: 






,,1 

(C4) 
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FIG. 1: Kinetic energy during the whole switch off of the hehum atom embedded in a 63 hydrogen 
system, in the 3-point quadrature case. The average kinetic energies are displayed in solid line, the 
instantaneous ones into dashed line. The vertical lines separate the different domains of constant 
switching parameter A. 
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FIG. 2: Total (solid line) and potential (dotted line) energies corresponding to the simulation of 
{63 H, 1 He} with the switching parameter equal to 0.5. 
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FIG. 3: Different values of ( rp,^\ ,^ a\ / obtained for the different quadratures. A quadratic fit 
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of the results is given as a guide for the eye. 
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